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1. Introduction 

Phenomenological laws of transport have been established long ago. In 1807 Joseph 
Fourier submitted his manuscript "Theorie de la Propagation de la Chaleur dans les 
Solides" to the French Academy. In it he presented several groundbreaking advances. 
He described heat conduction by a partial differential equation, equating time deriva- 
tive of the temperature with its laplacian, whereas previous attempts mostly focused on 
describing conduction between two discrete objects. By using a continuum description 
Fourier could use a linear equation (instead of some complicated "n-body"-like equation) 
for which he used a superposition principle. He solved the equation by separating vari- 
ables and using a trigonometric series for the solution - what is today known as Fourier's 
series. Fourier's law, stating that the heat current is proportional to the temperature gra- 
dient, lead subsequently to other phenomenological transport laws, like Ohm's or Ficks's 
law, all having the same form: current is proportional to the gradient of the driving field. 
According to standard Academy's practice, a referee committee consisting of Laplace, 
Lagrange, Monge, and Lacroix has been appointed in order to judge the suitability of 
the Fourier's manuscript. Laplace and Lagrange did not receive the manuscript very well 
though. In particular, they were not in favor of the trigonometric series that Fourier used 
in his solution and the manuscript has never been published by the Academy. Eventually, 
in 1822 Fourier published it by himself as "Theorie Analytique de la Chaleur". For more 
on rich history behind the Fourier's law see |[ll|2l- 

More than 200 years later microscopic origin of Fourier's law is still not very well 
understood |j3i]. Given a systems, for instance described by a Hamiltonian, one can ask 
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whether this system obeys Fourier's law or not. It is by no means granted that the Fourier's 
law will hold. For instance, taking a system of coupled harmonic oscillators Fourier's law 
is obviously violated. This happens because there are normal modes - the phonons - 
which do not interact and propagate in a ballistic fashion. In such a system the current is 
proportional to the temperature difference and not its gradient. Therefore, at fixed driving 
the current is in a ballistic system independent of the system's size n, j ^ nP, whereas 
it scales as j ^ 1/n if Fourier's law holds. Grossly speaking, the question therefore is, 
how does the current scale with the system size? If it is independent of n we say that 
the system is a ballistic conductor, if it decreases as '--^ 1/n we say that the system is 
a normal conductor (sometimes, by a small abuse of language, also diffusive, because 
diffusion implies Fourier's law; note though that Fourier's law does not necessarily imply 
diffusion). 

Based on the above example of a harmonic oscillator chain one could argue that all 
integrable systems would be ballistic. However, thing are not that simple. For instance, 
taking one of the simplest quantum models, a one-dimensional Heisenberg model, 

spin conduction for A > 1 is still debatable, see references in e. g. Heisenberg model 
has been introduced by W. Heisenberg in 1928 Isj] as a model system to explain the ferro- 
magnetism. Existing theories at the time, based on a direct interaction between magnetic 
moments, could not explain very high transition temperatures in ferromagnets. Heisen- 
berg used the symmetry of wave functions to derive the so-called exchange interaction, 
resulting in an effective Heisenberg model. Heisenberg model in Id is integrable by Bethe 
ansatz |6]. Despite integrability the evaluation of observables is rather involved. Calcu- 
lating time-dependent correlation functions, which could be used in the Kubo formula for 
conductivity, is at present not possible. Even for such simple system one therefore has to 
use either numerical computations of approximate techniques. If one uses Jordan-Wigner 
transformation |7] between spin-1/2 particles and spinless fermions one can rewrite Id 
Heisenberg model as a system of interacting fermions. It is therefore one of the simples 
so-called strongly-correlated systems, where many-body effects play a role. Although 
Heisenberg model has been initially thought of as being a simple toy model, today we 
know that it is realized to a very high degree of accuracy in many spin-chain materials, 
for a review see e.g |8]. 

The simplest transport situation is that of a stationary state, in which expectations do not 
change in time. Besides its relevance for transport such a stationary nonequilibrium state 
is one of the simplest nonequilibrium situations that one can study. Physics of nonequi- 
librium systems is much less understood that that of equilibrium. Detailed balance, which 
greatly simplifies equilibrium discussion, does not hold out of equilibrium. In stationary 
nonequilibrium situation only the total probability flow out of a state has to be zero and 
not along each connection individually as in equilibrium. Furthermore, there are very 
few exactly solvable nonequilibrium models from which we could learn about generic 
properties of nonequilibrium steady states (NESS). This is particularly true for quantum 
systems, whereas in classical physics there is a notable class of exactly solvable nonequi- 
librium stochastic lattice gasses |9]. In quantum domain there are basically to ways in 
which one can study stationary nonequilibrium phenomena: (i) one can consider Hamil- 
tonian evolution of a combined central-system + environment system and from its exact 
solution deduce the dynamics of the central system, or (ii) one can trace out environmen- 
tal degrees of freedom before solving the system, deriving a master equation describing 
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the evolution of the central system only. To use the first approach one obviously has to 
know how to exactly solve the total system. This is possible only in few instances, for in- 
stance for an XY spin chain 1 10] which is equivalent to free fermions, or for the so-called 
star system The second approach, in terms of a master equation, also allows for an 
explicit solution if the superoperator can be written as a quadratic function of fermionic 
operators 1 12]. 

In the present work we shall discuss quantum transport, and nonequilibrium physics 
in general, by studying solutions of master equation in a nonequilibrium setting. We are 
going to be in particular interested in the NESS as they are the simplest nonequiUbrium 
states. Master equation that we shall use is the Lindblad equation lfl4lfl5ll . 

^p = i[p,i/] + /:^-(p) (2) 

where the dissipative linear operator C^^^ is expressed in terms of Lindblad operators Lfc, 

fc 

Derivation of the Lindblad equation from the first principles, that is from a particular 
system-bath Hamiltonian, involves a number of approximations fi3\. Note though that 
Lindblad equation describes the most general completely positive trace-preserving dy- 
namical semigroup. Notwithstanding this, we know that there are evolutions that can not 
be described by a Lindblad equation. Our approach here to these issues is pragmatical. 
We shall not discuss approximations involved, neither if they are fulfilled, nor if Lindblad 
description is valid. The "philosophy" is that rather than study some more exact, but also 
more complicated, description about whose solutions we are not able to say much, we 
shall choose simple Lindblad equation with local dissipative terms which is amenable to 
solution. Our goal is to find solutions, even if they are for a simplified model system. 
After all, one of the reasons for the success of physics is that often simple model systems 
display generic behavior. In view of the lack of knowledge about quantum nonequilib- 
rium system we fell that description with local Lindblad equation is the way to proceed at 
this early state of investigations. 

We shall use two approaches to find the NESS of Lindblad equation. In the second part 
we shall present a model for which an exact solution is possible for any system size n. 
For systems which are not exactly solvable one has to either use approximations 1 16] or 
resort to numerics. Due to exponential scaling of computational complexity with n, until 
very recently, numerics has been limited to chains of fairly small size of less than « 20 



spins lllTl 11811 from which it is hard to conclude about properties in the thermodynamic 
limit which is of especial interest to us. Recently however, a time dependent density ma- 
trix renormalization group method (tDMRG) has been presented |19], which allows to 
study much larger systems. tDMRG can be used to simulate unitary pure state evolution 
as well as evolution described by master equation. Transient nonequilibrium studies of 
non-stationary pure states using tDMRG are quite abundant, an incomplete list includes 
for instance studies of transport via the spreading of inhomogeneities [20J or I — V char- 
acteristics ll2l[l . On the other hand there are much fewer studies of time evolution within 
a master equation approach. In our view master equation formulation has certain theo- 
retical and practical advantages over pure state evolution, see discussion below. In the 
first part we shall therefore present tDMRG method and its application on finding NESS 
and calculating transport properties. Both methods, analytical solution and numerical 
investigation, enabled us to study some interesting phenomena which are not found in 
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Figure 1. Schematic representation of our Lindblad equation l2ll. In addition to unitary 
part, dissipative part due to baths acts on the first and the last spin. 
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Figure 2. Schematic representation of a MPO with 7 sites. Each site is described 
by matrices (balls) having a local basis-label (Sj - yellow vertical lines) and two matrix 
indices (horizontal lines). 



equilibrium systems, like nonequilibrium phase transitions and a ubiquitous presence of 
long-range order. 

2. Numerical study 

We shall explain the working of tDMRG algorithm to calculate NESS on a case study 
of spin transport in the Heisenberg model ([T]i with the anisotropy A = 1.5 (J = 1) at 
infinite temperature reported in Il22ll . First, the dissipative Lindblad part C'^'^^ will act 
only on the 1st and the last spin. Each of these two parts has two Lindblad operators, 
one being proportional to (t+ = [a^ + icr^)/2 and the other to = {a^ — \(j'^)/2. 
If the coefficients in front of these two operators are different, they will try to induce a 
net nonzero magnetization on the corresponding spin. They can be expressed in terms of 



single bath parameter /i playing the role of driving potential Il22[l (in fact, /i = /?(/), if fi is 
the inverse temperature and 4> chemical potential). Having different ^ at the left and right 
ends, in our case we take /iicft = 0.02 and /bright = —0.02, will induce a nonequilibrium 
driving. Starting with some initial density matrix p(0) we want to calculate p{t) at time t 
being a solution of the Lindblad equation (|2|i. The asymptotic stationary state p{t — > 00) 
is then the sought-for NESS. To calculate p{t) as efficiently as possible we write it in 
terms of the so-called matrix product operator (MPO) representation, 

P = ^ E(l|^^'4"^ • ■ • A\:-\l) a^a^ (4) 

For each site j we have four matrices each of dimension D x D, for schematic 

representation see fig.|2] The advantage of a MPO representation lies in the fact that if the 
state p has small entanglement then the matrix size D can be small and independent of n 
so that the number of parameters needed to describe the state scales as ^ nD^, whereas 
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it would scale exponentially in n in the worst case of large entanglement because the 
matrices would have to be exponentially large. Interesting fact is that some states, like 
ground states in Id, do have such small entanglement. For more on matrix product states 
see the review [23]. 




Figure 3. Schematic working of one short time-step transformation exp (jCjj+iSt) 
acting on j-th and {j + l)-th spins of an MPO ansatz. Algorithm runs from top to 
bottom. 

Time evolution according to the Lindblad master equation now proceeds in small time 
steps of length St. The formal solution p{t) — exp {Ct)p{0) is decomposed using Trotter- 
Suzuki formula into small time-step operations on nearest neighbor spins. Decomposi- 
tion into operations acting only on nearest-neighbor spins is crucial and connected to the 
fact that matrix indices connect only neighboring spins, see fig. |2] Performing nearest- 
neighbor transformation on the MPO ansatz transforms the two neighboring matrices into 
one larger one. To restore a MPO form a singular value decomposition is performed, 
rewriting the state p{t + St) in a MPO form, but with the two matrices being of a larger 
size. If they have been of size D before the application of the transformation, they are 
now of size A - D. To prevent an exponential increase of matrix size with time we have to 
truncate them, usually back to the original dimension D. Doing that we of course make 
some truncation error The size of it depends on the discarded singular values, which 
are in turn equal to the squares of Schmidt coefficients, i.e. the entanglement. One short 
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time step of tDMRG algorithm, acting on nearest neighbors, is pictorially described in 
fig-E] Technical details can be found in original references 1 19], whereas our implemen- 
tation for the Lindblad equation is described in |22] and |24]. Crucial question for the 
efficiency of the algorithm now is how big must the matrices be for the truncation error to 
be bearable? For unitary time evolution (that is without a dissipative part) entanglement 
for generic (chaotic) systems always grows with time and with it also the necessary di- 
mension D. tDMRG method for conservative systems is therefore not efficient because 
it breaks down after relatively short time |25]. It seems that things are not that bad for 
open systems though. As a rule dissipation namely decreases the entanglement. Although 
no systematic study has been performed to date it seems that often one can sufficiently 
accurately describe NESS with matrices of the order D ~ 100 for chain lengths of order 
n ^ 100. Using tDMRG on the Lindblad master equation one can therefore calculate 
NESS for significantly larger systems than with other methods. 

For our spin conduction example we could calculate spin current jk = 2{a^aJ,_^^ — 
o'kO'k+i) '■^^ NESS for systems of various sizes. From the figure|4]one can see that the 
current indeed scales as jk ~ 1/n. At infinite temperature gaped Heisenberg model there- 
fore seems to display normal spin transport. Defining transport coefficient k in terms of a 
gradient of magnetization, j x K{zn — zi)/n, we get k k, 2.3. Approximately the same 
value of the transport coefficient is obtained also from the equilibrium correlation func- 
tion using Kubo formula Ii26il and qualitatively agrees with findings in 1 1 Sl 12711 . Because 
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Figure 4. Normal spin transport in Heisenberg model at A = 1.5; we show a scaled 
current divided by the magnetization difference. Horizontal line shows the value of the 
transport coefficient k ~ 2.3. Data from L221 . 



we have defined k in terms of magnetization, whose conjugate variable is /i = /30, k is 
actually the diffusion coefficient. Spin conductivity a, which is defined in terms of the 
gradient in the potential 0, is then simply a = /3k. 

For the end let us list some advantages and disadvantages of using tDMRG to calculate 
NESS. Among advantages are: (i) dissipative part of Lindblad equation seems to decrease 
entanglement significantly as compared to unitary evolution. Precise scaling is not known 
but simulating systems of rt « 100 spins is feasible with modest computational resources; 
(ii) For the efficiency of calculating NESS the convergence time after which we get the 
asymptotic state pit — cxd) is also relevant. The convergence time is given by the in- 
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verse gap of the superoperator £ (|2]l. The scaling of the gap has been studied in some 
exactly solvable systems and has been found to be always polynomial in n; (iii) 

Solving for NESS of a master equation one is actually as close as possible to an actual 
experimental situation. On the down side we can mention: (i) For stronger driving /i the 
entanglement increases and the method is not so efficient. This in particular makes it 
harder to study 12811 some interesting phenomena that occur only at strong driving, like 
negative differential conductivity, (ii) For density matrices tDMRG is more computation- 
ally intensive than for pure states (provided we have the same D). (iii) It works efficiently 
only in Id. (iv) No theorems exist about the entanglement of NESS states that would 
guarantee the efficiency of tDMRG as opposed to ground-state situations, (v) With local 
Lindblad operators it seems that it is in general difficult to simulate baths at very low tem- 
peratures 1 29]. In this regime it looks that pure state evolution using tDMRG 1,20,1 might 
have certain advantages. 



3. An exact solution 

An exact solution for NESS is possible in systems whose superoperator £ is a quadratic 
function of fermionic operators 1 12]. Recently an exact solution il3\\ has been found also 
for a system that is not quadratic in fermions. The unitary bulk part is described by XX 
Hamiltonian, 

Coupling to the baths at boundary two spins is the same is in the example in previous 
section. At each side we have two Lindblad operators. 



at the left end, while at the right end we have 



Lf = ^il + ti)a+, L^ = Vil^^n- (7) 

With only these terms the model is actually quadratic in fermions and has been solved 
before i30ll . An interesting twist now comes if we add an additional dissipative part due 
to dephasing at each spin. That is, each spin site experiences an independent reservoir 
that causes the decay of off-diagonal matrix elements. Such dephasing can be described 
by a single Lindblad operator at each site, 

i^r^=/|-|. (8) 

Dephasing makes the system nonquadratic because the dephasing part is quadratic in a^, 
which is itself quadratic in fermionic operators. This renders the system more interesting, 
for instance, the NESS is not gaussian and the Wick theorem does not apply lil3]. Nev- 
ertheless, the NESS can be calculated exactly in a systematic way ifisll . First, one notices 
that the equiUbrium stationary state, that is NESS for /i = 0, is a trivial identity, 

Peq = ^l- (9) 
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This can be interpreted as an infinite temperature state. For small driving /i the NESS will 
differ from pcq by a small amount. The ansatz for the NESS is written as a formal series 
in /z, 



-If 

2" V 



1 + Aii?^^^ + + . . . ^r^(r) + , . . ^ ^ (IQ) 



however, there is a crucial difference from the ordinary perturbation series. Due to al- 
gebraic structure, if one writes stationary equation Cp — Q satisfied by NESS, one can 
see that the unknown terms in the p-th order term i?'^' depend only on coefficients in the 
lower order terms. Now because the zeroth order term is explicitly known, it is 1/2", 
one can solve for terms order by order. Also, p-th order term is a sum of products of p 
factors, each one being either a magnetization a^, or a spin current jk- This also means 
that connected p-point correlation functions are contained solely in R'-p^ . First two orders 
are actually quite simple and we can write them in full. One gets 

n ^ n— 1 

fiR^^'> ^fiA + fiB, fiA^Y.^'i'^j' f^B=-Y,Jk- (11) 

j=i k=i 

Magnetization profile given by aj is linear and equal to 

2 + 2(j-l)7 + <5jM-J,-„ fi 

aj^~H + fi — — b=-—- — . (12) 

2+(n-l)7 2+(n-l)7 

We can see that the current scales as 1/n for large n. The system is therefore a normal 
conductor as long as the dephasing strength 7 is nonzero. 2nd order term can be written 

as 

„2 

^2^(2) ^ ^ ^y^^ ^ ^2(j ^ ^2jj ^ ^2p^ (.j3^ 



j=i k=]+i 




n+l-j 



f. 71 — 1 

kA = l 

For large n the dominant term in the 2nd order comes from Cj^k whose exact expression 
is 

Ci.fc = -7 o^ I 0/ + {n-l + 2/7)(5j+i,fc) , (15) 

[n- Z) + 2/7 

where Uj = 2 + 2(j — 1)7. Cj,k is equal to two-point connected correlation of magnetiza- 
tion. One can see that for 7^0 there is a long-range order, i.e., at a fixed distance \ j — k\ 
the correlation function goes to a constant plateau value of fi^/n independent of |j — k\. 
The plateau is of purely nonequilibrium origin and decreases in the limit of n 00. 
Similar long-range plateaus have been observed in quadratic nonequilibrium quantum 
systems 1 12], see also LSlil . leading to a conjecture li3Zl that this a generic situation. 
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4. Discussion 



Solvable model presented in the last section is in a way minimal model that can display 
normal transport. Any nonequilibrium steady state must necessarily have a nonzero cur- 
rent and magnetization terms, plus any optional correlations between these two quantities 
(these correlations are perhaps necessary to have nonzero gradient). In our solution these 
are exactly the terms present. Namely, the NESS is composed only of products of local 
magnetizations and currents. The model exhibits some other interesting properties. The 
NESS can be written in terms of a MPO ansatz with matrices of small dimension D = A 
that is independent of n il3il . It also displays a nonequilibrium phase transition at 7 = 0, 
going from a state without long-range correlations to the one with. 

On the numerical side, tDMRG shows to be a very useful method to explore nonequi- 
librium physics, especially in a stationary setting, where things are believed to be simpler 
than for instance in transient nonequilibrium situation (quantum quenches, for example). 
Considering that the field of nonequilibrium quantum physics is still rather unexplored it 
seems certain that many interesting things are still to be discovered. 
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